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Abstract. A new adaptive weighted essentially non-oscillatory WENO-# scheme in the con¬ 
text of finite difference is proposed. Depending on the smoothness of the large stencil used in 
the reconstruction of the numerical flux, a parameter 8 is set adaptively to switch the scheme 
between a 5th-order upwind and 6th-order central discretization. A new indicator r 6 measuring 
the smoothness of the large stencil is chosen among two candidates which are devised based on 
the possible highest-order variations of the reconstruction polynomials in L 2 sense. In addition, 
a new set of smoothness indicators s of the sub-stencils is introduced. These are constructed 
in a central sense with respect to the Taylor expansions around the point Xj. 

Numerical results show that the new scheme combines good properties of both 5th-order up¬ 
wind schemes, e.g., WENO-JS ( JJS96I |1. WENO-Z f |BCCD08] 1. and 6th-order central schemes, 
e.g., WENO-NW6 ( IYC09I 1. WENO-CU6 f lHWAlOl b In particular, the new scheme captures 
discontinuities and resolves small-scaled structures much better than the 5th-order schemes; 
overcomes the loss of accuracy near some critical regions and is able to maintain symmetry 
which are drawbacks detected in the 6th-order ones. 
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1. Introduction 

In this work, we consider the following one-dimensional hyperbolic conservation law 


( 1 . 1 ) 


where u = (u\, 


u t + f(u) x = 0, x £ 
u(x, 0) = U 0 (x), 


t > 0, 


, u m ) T is an m-dimensional vector of conserved quantities and its flux f (u) is 
a vector-valued function with m components, x and t denote space and time, respectively. Eq. 
fll.lH is called hyperbolic assuming that all eigenvalues A s ’s of the Jacobian A(u) = df/du are 
real and the set of all eigenvectors r s ’s is complete. 


It is well-known that shocks and discontinuities may develop in the solution of Eq. (1.1) even 


if the initial condition is smooth. Thus classical numerical methods which depend on Taylor 
expansions in general do not work in this case. As a result, there exist spurious oscillations near 
these discontinuities. 
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In order to overcome this difficulty, in |Ha83] . (Ha84) Harten introduced the Total Variation 
Diminishing (TVD) schemes which are of high-order resolutions as well as oscillations free. The 
schemes are constructed based on the principle that the total variation of the numerical approx¬ 
imation must be non-increasing in time. A drawback is that TVD schemes are only at most 
first-order near smooth extrema (see [OC84] L Later on, Harten et al. in |HQEC86j . |HQ87j . and 
(HEOC97| tried to tackle this disadvantage by relaxing the TVD condition and allowing spuri¬ 
ous oscillations in the order of the truncation error to occur but the 0(1) Gibbs-like ones are 
essentially prevented. Thus these new schemes were named essentially non-oscillatory (ENO). 
For an rth-order ENO scheme, only the smoothest stencil is chosen among r candidates to 
approximate the numerical flux. The smoothness of the solution on each stencil is determined 
by an indicator of smoothness. Later on, Liu, Osher, and Chan ( jLOC94j ) upgraded ENO 
schemes and introduced the Weighted ENO (WENO) by combining all stencil candidates (here¬ 
after sub-stencils) in the numerical flux approximation. Here, a nonlinear weight is assigned 
to each sub-stencil to control its contribution in the procedure. WENO schemes maintain the 
essentially non-oscillatory property of the ENO near discontinuities and outperform the latter 
in smooth regions where the accuracy order is increased to (r + l)th-order if r sub-stencils 
are used. Consequently, Jiang and Shu (see (JS96 ]. also [Sh03], and the review [Sh09| l con¬ 
structed WENO schemes in the framework of finite difference and further improved the order 
to (2 r — l)th in smooth regions by introducing a new class of smoothness indicators. Hereafter, 
we denote WENO-JS for the 5th-order finite difference WENO developed in [JS96I . In [ BSOO j. 
j'SZ'()~8| higher order than 5th-order WENO schemes are given. 

Since the introduction of WENO, many improvements and derivatives of the schemes have 
been developed and introduced. Henrick et al. in [HAP05] carefully analyzed the necessary and 
sufficient conditions of the nonlinear weights and found that WENO-JS does not achieve the 
designed 5th-order but reduces to only 3rd-order in cases where the first and third derivatives 
of the flux do not simultaneously vanish (e.g., f'(xj) = 0 but f"(xj) ^ 0 for the scalar case 
of Eq. &■ They then suggested an improved version which is called mapped WENO, 
abbreviated by WENO-M. By using a mapping on the nonlinear weights, WENO-M satisfies 
the sufficient condition on which WENO-JS fails and obtains optimal order near simple smooth 
extrema. In a different approach on the construction of the nonlinear weights, in |BCCD08] 
Borges et al. introduced the 5th-order WENO-Z scheme. Here, the authors also measured 
the smoothness of the large stencil which comprises all sub-stencils and incorporated this in 
devising the new smoothness indicators and nonlinear weights. It was proven numerically that 
WENO-Z is less dissipative than WENO-JS and more efficient than WENO-M, respectively. 
It was also checked that WENO-Z attains 4th-order near simple smooth extrema comparing 
with 3rd-order of WENO-JS. For higher order WENO-Z schemes, we refer readers to [CCDllj . 
Another approach to improve WENO schemes is the new designs of the smoothness indicators. 
In |HKLY13| . L 1 -norm based smoothness indicators are suggested, and the ones devised from 
Lagrange interpolation polynomials are given in |Fal4j , and jFSTY14] . See also jFHW12] for a 
new mapped WENO scheme. 

We notice that for a general flux where the signs of the eigenvalues of the Jacobian A(u) are 
not uniform throughout the domain, a flux splitting technique, for example, the global or local 
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Lax-Friedrichs or the Roe with entropy fix (see }JS96j and the references therein) is needed. 
This increases the number of grid points in the numerical flux approximating procedure by one. 
We take the 5th-order WENO-JS scheme for example, the total number of grid points used 
in the reconstruction for both positive and negative fluxes will be six instead of five. We also 
note that with these six points, one can indeed improve the scheme up to 6 th-order in smooth 
regions. The difficulty of this approach lies in the dispersive nature of a central scheme if six 
points are employed. In this case, oscillations are expected to occur near discontinuities. In 
(see also |CFY13] for the boundary condition treatment), Yamaleev and Carpenter for 
the first time introduced a 6 th-order WENO scheme by adding one more sub-stencil into the 
numerical flux approximation. We denote this scheme WENO-NW 6 . For this most downwind 
sub-stencil, an ad hoc treatment on the smoothness indicator was suggested. The idea is 
originated from that of Martin et al. in jMTWOfij . In order that oscillations do not happen, /I 3 
is computed using the information of f(u) on all grid points of the large stencil, i.e, six points. 
Hence, the sub-stencil only plays roles in case the solution is smooth over this large stencil. In 
a similar manner, recently Hu, Wang, and Adams in IHWAlO] proposed an adaptive central- 
upwind WENO-CU 6 scheme which switches between a 5th-order upwind and 6 th-order central 
WENO scheme automatically. The difference of their work from that given in |YC09] is that 
@3 is defined via a Lagrange interpolating polynomial of degree five over the large stencil. In 
|HA11| , the authors successfully applied WENO-CU 6 in the LES simulation of scale separation. 
Other hybrid WENO schemes can be found in, for examples, jCD07j . or |LQ10| , |HP04l . etc. 

A drawback of the presented 6 th-order WENO schemes (i.e., WENO-NW 6 , WENO-CU 6 ) is 
that they suffer from a loss of accuracy near the smooth critical region which is just behind 
another one where the first derivative of the flux is undefined. To illustrate this, we consider 


Eq. (1.1) in a scalar case where f(u) = u in the following example. 


Example 1.1. 


( 1 . 2 ) 


u t + u x = 0, x (E (-1,1), 
uo(x) = max(—sin(7rx), 0), 
subject to periodic boundary conditions. 


We approximate the solution of (1.2) by the WENO-JS, WENO-Z, WENO-NW 6 , and WENO- 


CU 6 schemes. The results at time t = 2.4 with 200 grid intervals are plotted in Fig. [l] with 
the critical region zoomed in. It is clearly shown the above mentioned defect of the WENO- 
NW 6 and WENO-CU 6 schemes. Near the smooth critical region, we note that these schemes 
are worse than both WENO-JS and WENO-Z. Since there are many problems whose solution 
often exhibits the same behavior as mentioned above, we notice that this loss of accuracy is an 
important issue. 


Our goal in this work is to construct a new WENO scheme which overcomes the drawback 
of WENO-NW 6 and WENO-CU 6 presented in the previous example. For this, we introduce 
a different switching mechanism between a 5th-order upwind and 6 th-order central scheme. 
Unlike the WENO-NW 6 or WENO-CU 6 scheme in which the change depends on the smoothness 
indicator of the most downwind sub-stencil, in our scheme, whether the scheme is upwind or 
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FIGURE 1. Left: Numerical solutions of Eq. (1.21 at time t = 2.4 obtained from 
different WENO schemes. Right: Zoom near the critical region. 


central is due to the smoothness indicator of the large stencil. Moreover, instead of using all 
six points for the indicator fy, we reduce the number of points down to only four. The reason 
for this is explained in the below section. We also introduce a new set of smoothness indicators 
which are constructed in a central sense in Taylor expansions with respect to the point Xj. The 
feature of our new scheme is that it automatically switches between a 5th-order upwind scheme 
near discontinuities to prevent spurious oscillations, and a 6th-order central scheme in smooth 
regions which improves the loss of accuracy of the WENO-NW6 and WENO-CU6 schemes. 
Moreover, it is shown in the numerical results below that the new scheme maintains symmetry 
in the solutions much better than the 6th-order ones. 

We organize our paper as follows. In section 2, we summarize the mentioned above WENO 
schemes which relate to our work. From there, we construct our new scheme in section 3. In 
this section, we first start with the new definition of the central smoothness indicators. We 
then introduce a new switching mechanism for 5th-order upwind and 6th-order central scheme. 
Numerical results comparing the performances of all schemes are presented in section 4. Finally, 
we close our discussions with a conclusion section. 


2. Summary on Finite Difference WENO Schemes 


For simplicity, we consider Eq. (1.1) in a scalar case. We rewrite the equation as follows, 


( 2 . 1 ) 


u t + f(u) x = 0, x E [xi, x r \, 
u(x, 0) = Uq(x). 


We first define a uniformly spatial grid Xj = x;+jAx, j = 0,. .., N, where Ax is the grid size. 
We denote the interval Ij = [x ■_ i, x - + i] where x- ± i = Xj ± are the interfaces of Ij. We also 
denote all quantities with a subscript (-)j their grid values at Xj, for examples, Uj = u(xj,-), 
fj = etc.; and so as with a subscript (-) . + i for the quantities at the interface x- + i. 

Whether these quantities are exact or approximate depends on particular circumstances. 
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( 2 . 2 ) 


We denote h(x) the numerical flux function defined as follows, 

. A n ? 

1 


f(u{x, •)) = h(x) ■= -z- ( h{y)dy. 

Ax J x _^l 


Evaluating Eq. (2.1) at grid point Xj. we obtain the semi-discretized form as follows, 

rfoi . f) f 

(2.3) 


dt 


3/ 

dx 


K x j+i) - 


Ax 


=: C(u). 


It is noticed that Eq. (2.3) is exact since there are no approximating errors in the formula. 


2.1. Time Integration. We first mention about time advancing for Eq. (2.3). Following 


|GS98] and the references therein, for all below WENO schemes, we employ the 3rd-order TVD 
Runge-Kutta method as below. For TVD, we mean that the time integrator follows the same 
property mentioned above. Here, At is the time step satisfying some proper CFL condition. 

- (1) = u n + AtC(u n ), 


(2.4) 


u K 


u (2) = -y> + V 1 ) + -a *£(«W), 

4 4 4 v ' 


= -u n + -u'. 
3 3 


,n+l _ ± n ,n ^ ±„(2) + _AtT( n ( 2 )), 


where C{u) obtained from some method is an approximation of the spatial operator C(u). In 


•, i follows 

J 2 


q. (2.47) with q^’s replaced by 


particular, see the below WENO discretizations £ 5 (u) in Eq. (2.12) where f 5 i = f, 

_ ^ ^ __ J+2 

Eq. (2.13) and C e (u) in Eq. (2.46) where / 6 , is defined in 
_ _ -7+2 

Wfc’s. 

We now proceed to the discussions on the spatial discretizations. 

2.2. 5th-order Upwind WENO Reconstruction. We notice that for simplicity, we can 
assume that f{u) > 0 over the whole computational domain. In case there is a change in signs 
of f(u), a flux splitting technique is invoked. We discuss this in Remark 


2.1 


below. 


Originally, WENO schemes were constructed in the context of finite volume (see jLOC94j h 
Thanks to Lemma 3.1 given in |Sh09] . the schemes can be transformed into finite difference 


through relation (2.2). The hj is called an average value of the numerical flux h{x) over the 
interval Ij. We then seek for an approximating polynomial y(x) of degree four of h(x) as below 

(2.5) h(x) ~ f 5 (x ) = ao + a\x + a^x 2 + a^x 3 + a^x A , 

over the large stencil S 5 = {xj- 2 , Xj-%, Xj, Xj + 1 , Xj+ 2 }- We note that S' 5 is chosen biased to the 

left with respect to the point x -, 1 for the stability purpose. Hence, the scheme is in an upwind 

■? i_ 2 

sense. 


Replacing the integrand h{x) in Eq. (2.2) by its approximation f 5 (x) in Eq. (2.5) and 
evaluating at xk = j— 2,... ,j+ 2, we can uniquely determine the coefficients s, k = 0,..., 4. 


Since the procedure is via the average hk = h{xk) = f(u(xk,-)) in Eq. (2.2), it is called 
reconstruction and f 5 (x) is the reconstruction polynomial. Evaluating f 5 (x) at x- + i , we obtain 
the approximation of h ■, 1 as follows, 

J' l_ 2 

£5 = l f . _ 13, ,47 27 

U+5 60 /j 2 60 /j 1 + 60 /j + 60 /j+1 60 Ij+2 ' 


( 2 . 6 ) 
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X j-1/2 


X j+1/2 



X. . 

X. , 


X. 


X. , 

X. _ 
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H 

J + 2 

□- 



—□ 





S 


1 


Q- 


■B 


S 


2 


Q- 


■B 


FIGURE 2. Stencils for 5th-order WENO schemes. 


Here, we recall f k = f(u k ) = f(u(x k , ■))• 

To justify the approximation error, we denote the polynomial 77(x) such that 


(2.7) 


We then deduce from Eq. \2.2\ that 

( 2 . 8 ) 


f(u(x,-)) = 


1 


x+^L 


H'{x) = h(x). 

H{x + %E) - 77(x - ^ 


Ax 


H'(y)dy = 


Ax 


Substituting Eq. (2.8) into the approximation (2.6), evaluating at x k , k = j — 2,..., j + 2, 
and applying Taylor expansions of 77(x) at x = x-,i, we obtain the following truncation error, 

J-|- 2 


(2.9) 


7h = 


. (—277 -5 + 1577-3 — 60i7_i + 2077 ,i + 3077 ,3 — 377 ,5 ) 

oOAx ^2 J 2 J 2 ^ ' 2 ^2 2 


H , 1 d 6 77 

i+l 60 dx 6 


Ax 5 + 0(Ax 6 ) = h 1 - 


1 d 5 f 


x=x 1 
■ 2+2 


J+ 2 60 cte 5 


Ax 5 + 0(Ax 6 ). 


The last equality in Eq. (2.9) is justified as follows. Thanks to the relation in (2.8), by a 
Taylor expansion around x ■, 1 we have that 

•7 ' 9 


+ C>(Ax). 


d 5 f 

1 

d 5 H 


d 5 H 


d 6 H 

dx 5 

Ax 

dx 5 

x=x 1 

J-t 2 

dx 5 

x=x . , 1 — Ax 
7+2 J 

dx 6 


Together with Eq. (2.7), we deduce the last equality in Eq. (2.9). 
Similarly, we have 


( 2 . 11 ) 

Hence, we have 

( 2 . 12 ) 


p = A f _13, ,47 27 3_ 

L-5 60 /j 3 60 /j 2+ 60 /j + 60 /j 60 /j+1 


= h A AZ 

7“ 2 60 <9x 5 


Ax 5 + C>(Ax b ). 


diii 


i 5 7 5 

■'7+3 ■'j-3 


dt Ax 

The scheme is 5th-order of accuracy in space. 


= : £ 5 (U). 
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For non-smooth solutions, we employ WENO reconstruction. The idea of WENO schemes is 
that, instead of the 5-point stencil 5 s , a convex combination of three 3-point sub-stencils are 
facilitated for an adaptive choice of candidates for the reconstruction. That is, 

2 

(2.13) 


I' 


k=0 


where f k k ’s are defined below and ujk is the non-linear weight satisfying ujk > 0, V/c and 
•f'2 


(2.14) 


y = i. 


fc =0 


The necessity of non-negative nonlinear weights is discussed in |LSZ09| and in [ SHS02 ] for 
practical implementations. And f k + i is the approximation of h- + i by the reconstruction poly¬ 
nomial f k (x) = bo + b\x + b 2 X 2 over the sub-stencil Sk, k = 0,1, 2. Here, Sq = {xj- 2 ,Xj-i,Xj}, 
S\ = {xj-i,Xj,Xj+i}, and S 2 = {xj,Xj + i,Xj + 2 } (see Fig. [ 2 ]). Carrying a similar process as for 
the large stencil S 5 , we find that around xj = 0, 


(2.15) 

(2.16) 
(2.17) 


f°(x) = 


-fj-2 + 2/j-i + 23 fj 


24 


f\x) = 


-fj -1 + 26 /j ~ fj +1 

24 


fk) = 


23 fj + 2 f j+ i - fj + £ 

24 


fj -2 ~ 4/j-i + 3/j 
2Ax 

fj± 1 ~ fj-i 

2Ax 

~3 fj + 4/j+i - fj +2 
2Aa; 


f )—2 - 2/j-l + /; 


2Aai 2 


J 'x 2 , 


fj-i - 2fj + /j+i A 2 
2Aa: 2 / ’ 


fj ~ 2 fj+i + /j +2 \ 2 

2As 2 ' ‘ 


Evaluating each of these f k (x)'s at x -, 1 , we obtain that 


(2.18) 

fk 

(2.19) 

f +i 

(2.20) 

ki 


J+3- 

2 „ 7 . 11 


6 " 


6 


>3- 


6 


Jji 


2 

6 ‘ 


5 

6 " 


1 

6 " 


that 

( 2 . 21 ) 


Carrying a similar process as given in Eqs. (2.9) - (2.10) with f 1 replacing f 1 , we obtain 

«? ' 2 


V|=^ + °(A- 3 )- 


the following linear relation 
( 2 . 22 ) 


Comparing between f , given in Eq. (2.6) and the ones in Eqs. (2.18) - (2.20), we deduce 

f> + 2 


fj+ 


k=0 


where 


1 


6 


3 
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are called the linear (optimal) weights. We note that 


(2.24) 


= I- 


k =0 


Adding and subtracting Ylk= 01 'k.f ' .i into and from Eq. (2.13), thanks to the truncation 

_ _ •?' 2 1 _ 

errors in Eqs. (2.9), (2.21), the normalization in Eqs. (2.14), (2.24), and the linear relation 

( 2 . 22 ) we obtain that 

4ti = )/j±i 


k =0 


k =0 


?± i 
J 2 


(2.25) 


XX* -7fc)(V| +°( Ax3 )) + ( ~ 

fc=o “ V 


1 3 5 / 


J± 2 60 <9x 5 


60 dx 5 

Ax 5 + ^(cn± - 7 ±)0(Ax 3 ) + 0(Ax 6 ) 


Ax 5 + 0(Ax 6 


fc =0 


■d±§ 


where 7 ^, u;± are the linear and non-linear weights of the sub-stencils S J k " 2 corresponding to 

the interfaces x-ii, respectively. 

-' =c 2 _ 

Hence, in order that the discretization in Eq. (2.12) where / 3 ±1 = fj±i follow the nonlinear 
relation (2.13) to be 5th-order, we deduce the sufficient condition for the nonlinear weights as 
follows, 


(2.26) 


-7fc = 0(Ax 6 ), Vfe. 


Different WENO schemes depends on how these nonlinear weights and the smoothness in¬ 
dicators are defined. The latter ones are introduced in the below section. In the following 
subsections, we summarize the 5th-order upwind and 6 th-order central WENO schemes dis¬ 
cussed previously. Since WENO-Z is a good replacement for WENO-M, we omit the latter one 
in our comparison. 

2.2.1. WENO-JS. In |JS96] . Jiang and Shu defined the nonlinear weights as follows, 


(2.27) 


n J S 

JS a k 

k ^2 is ’ 


Xu=0 a i 


= 


Ik 


(£ + Pk) p ' 


where 7 *, is defined in Eq. (2.23), /3k is called the smoothness indicator of Sk which measures 
how smooth the solution is over this sub-stencil. The authors defined these /3k s through the 
normalized L 2 -norm of high-order variations of the reconstruction polynomials given in Eqs. 


(2.18) - (2.20). Explicitly, for a 5th-order scheme, we have 

1-2 


(2.28) 


/3 k = Ax 


A - 1 

3 


2 ( d 2 

\dx 2 


f k (x) dx, 


where / fc (x)’s are as in Eqs. (2.15) - (2.17) and the sub-stencils Sq, Si, S 2 are centered around 


Xj = 0 . 
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Evaluating for each k, we obtain that 


(2.29) 


A) - 2 _ 2 /?-i + fj ) 2 + ^(/j-2 ~ 4 /?—i + 3 /i) 2 

= / /2 Ax 2 + (^/" 2 - Ax 4 + 0(Ax 5 ), 


(2.30) 


IQ 1 

A = Y2 ^ i_1 ~ 2 ^' + ^ +1 ) 2 + 4 (/7+ 1 ~ fj-i ) 2 

= f 2 Ax 2 + (^r 2 + Ax 4 + 0(Ax 6 ), 


(2.31) 


A _ 12 ~ 2 h +l + ^+ 2 ) 2 + 4(^/7 - 4/j+i + fj+ 2 )' 

= f' 2 Ax 2 + f^f" 2 - p'f") Ax 4 + 0(Ax 5 ), 


where the derivatives are evaluated at x = Xj. 


In formula (2.27), e is a small parameter to prevent division by zero. In most cases, WENO- 
JS works well with e = 10 -6 . A thorough analysis of the role of e can be found in [HAP05) . 
The parameter p is to increase the dissipation of the scheme. For WENO-JS, p = 2 is chosen. 
If fj = f(xj) f 0, Vfc, /3k can be written in the form 


(2.32) 


Ac = (/, / Ax) 2 (l + 0(Ax 2 )). 


Substituting these into Eq. (2.27) with the removal of e, since (1 + y) 2 = 1 + 0(y), by Eq. 


(2.24) we obtain that 


(2.33) 


, ,JS _ 
u k — 


Ik(f'jAx) 2 (l + e>(Ax 2 )) 

(/;Ax)- 2 £Lo7K1 + 0(A* 2 )) 


= 7 fc + 0{ Ax 2 ), 


which is a relaxed form of (2.26). We notice that for WENO-JS, uji s cannot satisfy con¬ 


dition (2.26) directly. Moreover, if f- = 0, it is observed from Eqs. (2.29) - (2.31) that 


/3k = j|(/j') 2 Ax 4 (l + 0{Ax)) for k = 0,2 and (3\ = ^|(/") 2 Ax 4 (l + 0(Ax 2 )), in which the 


condition ( |2.32 ) is not satisfied for all AJs. Similarly, we find that 
(2.34) u)( s = 7 fc + O(Ax), 

which is a loss of accuracy near this critical point. This accuracy loss is improved by the 
WENO-Z scheme which is summarized in the next subsection. 


2.2.2. WENO-Z. Borges et al. in [BCCD08] proposed a new WENO-Z. In their scheme, the 
nonlinear weights are defined differently from those of WENO-JS. They are as follows, 


(2.35) 


= 


at 


Z ’ 

2^i=o a i 


= Ik 1 + 


T 


£ 3- (3 k 


where the smoothness indicators (3k s are the same as those given in Eqs. (2.29) - (2.31), 
e = 10~ 40 , and t z is the smoothness indicator of the large stencil S 5 . The power q is used to 
tune the relation between the dispersive and dissipative property of the scheme. It is checked 
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numerically in |BCCD08j that the scheme becomes more dissipative when q is increased. For 
WENO-Z, t z is defined as follows, 

13, 


= i A) - &i = f \rr\Ax 5 + o{ ax 6 ). 


(2.36) 

We note that if /j A 0 and q = 1, f3k = 0(Ax 2 ), \/k. Then 

(2.37) 


— = 0( Ax 3 ), Vk. 
Pk 


Similarly to Eq. (2.33), we obtain that 
(2.38) = 7fc + 0(Ax 3 ), 


directly without using the relaxed version as the WENO-JS scheme. It was also proven in 
[BCCD08j that WENO-Z is 4th-order near simple smooth critical points (i.e. where /j = 0) 
for <7 = 1 and attains the designed 5th-order for q = 2. The tradeoff for the latter case is that 
the scheme is more dissipative. Throughout this work, we choose q = 1. One more advantage 
of WENO-Z over WENO-JS is that the former is more central in a sense that the stencil over 
which the solution is discontinuous plays more roles in the approximation of the numerical flux. 
This assessment is checked as follows. We suppose that S 2 contains a discontinuity whereas 
the solution is smooth over the other two sub-stencils. Hence, t z = 0(1), /?2 = 0(1), and 
/3k = 0(Ax 2 ), k = 0,1. Then, 


(2.39) 


af 72 (l + g) 

“* T* (‘ + m) 


g(& + T Z ) = a J 2 s {p2 + T z ) of 
f k (pk + T z ) a( s (Pk + t z ) - a( s 


since 


(2.40) 


f3 2 + T Z ~ P 2 + T z 

j3k + T z T z 


> 1 . 


z — 

Hence, WENO-Z has a sharper capturing of discontinuities than WENO-JS. 


Remark 2.1. 

i. In case the condition f'{u ) > 0 is not satisfied, which is general in real applications, we 
apply a flux splitting technique to decompose f(u) into positive and negative compo¬ 
nents. In most applications, the global Lax-Friedrichs flux splitting is used (see [JS96j 
and the references therein), 

(2.41) / ± (^) = ± oiu), 

where a = max | f'(u)\ over the whole computational domain. Then, 


(2.42) 


/ 


J+5 


= f k +f i+i- 


The negative flux / . 

^+2 


is symmetric to f + , 
J+2 


'2 ’ 


with respect to x -, 1 

If the flux splitting is employed, the overall number of grid points used in the recon¬ 
struction of the numerical flux is increased by one. That is, let S' 5 " 1 " and S 5 ~ be the 


stencils over which f + , 
J+5 


and f. , are determined, then 


S 6 := S’ 5 + |j5 5 = {xj-2,Xj-l,Xj,Xj + l,Xj + 2,Xj +3 }, 


(2.43) 
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>2 
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J*3 


S 


1 


B- 


B 


5 2 B-B 

5 3 B-□- 


B 


FIGURE 3. Stencils for 6th-order WENO schemes. 


which consists of six points. We note that with this S' 6 , there exists a polynomial of 
degree five which reconstructs h{x ) over the stencil. Therefore, the accuracy of WENO 
schemes can be increased up to sixth. These schemes are discussed in the subsection 
below. 


2.3. 6th-order Central WENO Reconstruction. Carrying a similar procedure as described 
in the previous subsection for the 6 -point large stencil S 6 , we can deduce that 


(2.44) 

Similarly, 

(2.45) 


fj+b un ^ 7 2 1 60^ 60 ^ +1 60^ +2 60^ J+3 

„6 


60 
i d 6 / 
J+ 2 140 dx 6 


Ax 6 + 0(Ax 7 ). 


X=Xj 


37 


37 


4-1 60 ^ 3 60 ^ 2 + 60 ^' 1 + 60 ^' R0^ +1 + 60^ +2 


60 


. , J_3V 

140 cte 6 


60 2 60 J 60- 
Ax 6 +C(Ax 7 ). 


60- 


X = Xn 


Hence we obtain that 


(2.46) 


dUn 


76 76 

Jj+l 

3 2 3 2 =:C 6 {u). 


dt Ax 

The scheme is increased to 6 th-order of accuracy in space. 

Adding one more sub-stencil S 3 = Xj + 2 , x j+:i} into the approximation of the interfaced 

value fj + i (see Fig. 


3), we deduce a similar linear relation with Eq. (2.22) as follows 

3 


(2.47) 




k =0 


where 
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Here, / 3 1 is the 3rd-order approximation of the numerical flux h(x ) at the interface x 


3 +_ 

from the reconstruction polynomial / 3 (x) over the sub-stencil S 3 . Explicitly, we have 

11 . 7 „ 2 

6" 


j+i 


(2.49) 


fj -ft g /?+l ^/?'+ 2 + R fj+ 3- 


6 ‘ 


The other approximations f k 1 ’s, k = 0,1, 2 follow Eqs. (2.18) 

■? + 2 


( 2.20 . 


The nonlinear combination using the nonlinear weights is similar to that of the linear case 


(2.47), except for the linear weight 7 ^ replaced by its nonlinear version k = 0,..., 3. 


Remark 2.2. 

i. Since the accuracy order is increased to sixth, the sufficient condition for the nonlinear 


weights given in (2.26) is also increased by one. That is, 


(2.50) 


Wl 


- 7 , ± = 0(Ax 4 ), Vfc. 


11. 


Observing from Eq. (2.47) that the approximations /■_, i’s are now symmetric with 


respect to x-,i. It means that the scheme now is central. Hence, spurious oscillations 
are expected to occur near discontinuities. A treatment on the most downwind sub¬ 
stencil S 3 is needed to sustain the essentially non-oscillatory property of the scheme. 
We now overview the 6 th-order WENO schemes for this case. 


2.3.1. WENO-NW 6 . In IYC09I , Yamaleev and Carpenter proposed a 6 th-order energy-stable 
WENO scheme. They introduced an artificial dissipative term and proved that this makes 
the new scheme be stable in L 2 sense. In this work, we only discuss their treatment on the 
nonlinear weights and omit this artificial dissipative term (see |YC09] for a detailed discussion 
on the term). Hence, the scheme here is denoted by WENO-NW 6 , not ESWENO as in their 
paper. 

The nonlinear weights (denoted ) follow those defined in the WENO-Z scheme 


which are given in Eq. (2.35), for k = 0,..., 3. The differences lie on the smoothness indicator 


of the most downwind sub-stencil /I 3 and the one for the large stencil S 6 . For the former, in 
order that there are no oscillations occurring near discontinuities, all grid values of the flux over 
S 6 are accounted for the computation of (3 3 . It is as follows, 


(2.51) 


/?3 — |(/3q + fit + P 2 + ^l) 1//4 ) 


where /3% is computed using the formula given in Eq. (2.28). That is, 


(2.52) 


_ 12 ~ 2 ^'+ 2 + /?+3 ) 2 + ^( - 5 /j+ 1 + 8 fj +2 - 3/j+3)' 


= f ' 2 Ax 2 + ( ^/" 2 - ) Ax 4 + 0{ Ax 5 ), 


13 , a _ 11 
12 J 3 

where the derivatives are evaluated at x = Xj. The other indicators /Vs, k = 0,1, 2 follow Eqs. 


(2.29) - (2.31). 


The smoothness indicator of the large stencil S e is defined as the highest, i.e., fifth-degree, 
undivided difference as follows 

t nw = (/,-— 2 - 5/7-1 + IO/7 - 10/7+1 + 5/7+2 - fj + 3 ? 

= (/ ( 5 ) ) 2 Ax 10 + O(Ax n ). 


(2.53) 
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Carrying a similar procedure as in Eqs. (2.37) - (2.38) with a change in r, thanks to Eq. 


(2.53), we obtain that 


(2.54) 


r NW 

fik 


0(Ax 8 ), if/'/O, 

0(Ax 6 ), if fj = 0 , 


thus condition (2.50) is satisfied. Hence, the scheme is 6 th-order in smooth regions. The case 
where the derivatives of / vanish will be checked numerically in the below section. 


2.3.2. WENO-CU6. In |HWA10j . Hu et al. developed the adaptive central-upwind scheme 
WENO-CU 6 based on the principle that the most downwind sub-stencil only plays roles in 
smooth regions and is suppressed near discontinuities. Hence the scheme is central in smooth 
regions and upwind near discontinuities. The scheme is different from WENO-NW 6 in defining 
the smoothness indicators 03 and r. In particular, they are as below, 


(2.55) 

!h = ^^[271779/jU + /t-2(-2380800/,_i +4086352/,- - 3462252/, +1 + 1458762/, +2 - 245620/,- +3 ) 
+ /,_i(5653317/,■_! - 20427884/,- + 17905032/,- +1 - 7727988/, +2 + 1325006/, +3 ) + /,(19510972/, 

- 35817664/,+i + 15929912/,- +2 - 2792660/,- +3 ) + f j+1 (17195652 /, +1 - 15880404/,- +2 
+ 2863984/, +3 ) + /, +2 (3824847/, +2 - 1429976/, +3 ) + 139633/J +8 ] 

= f 2 Ax 2 + ^/" 2 Ax 4 + 0(Ax 6 ). 


It is noticed that there is a typo in Eq. (25) in IHWA10I . and we have corrected it in Eq. \2.55fy . 
From there, the smoothness indicator of the large stencil S 6 is defined as follows, 


(2.56) t cu = 0 3 - |(/3 0 + 40i + 0 2 ) = 0(Ax 6 ). 

Hence, we have for k = 0,..., 3,, 


(2.57) 



0( Ax 4 ), if/'/O, 

0(Ax 2 ), iff' =0, 


which satisfies the condition (2.50). 


It is also noteworthy that a )+ in Eq. (2.35) has a change as below, 

n-CU 


(2.58) 


a k U —7 k[C + 


£ + Pk 


where C 1 is to increase the contribution of the linear weights when the smoothness indicators have 
comparable magnitudes (see ITWM07] ). Following (HWA101 . we choose C = 20. 

As indicated in example |l.l[ the 6 th-order WENO-NW 6 and WENO-CU 6 schemes suffer from the loss 
of accuracy near the smooth critical points just right behind, with respect to the characteristic direction, 
a critical point where the first derivative of the solution is undefined, that is, the solution is just C° at 
that point. The explanation for this defect is given in the below section. In the next section, we propose 
a new scheme which automatically switches between a 6 th-order central and 5th-order upwind scheme 
and overcomes the defect occurred in the mentioned schemes. 
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3. The New Scheme 


We first observe that the 5th- and 6tli-order linear approximations given in Eqs. (2.61 and (2.44) 
respectively, can be combined linearly in the following manner, 


(3.1) 
where 

(3.2) 


/;+* =7o fl /? +i +T*/7 +i +71/^4. 


jj + 


i+3 


'1 + 3 


7 » = l (1 + ()) , 7? = 1(3 + 0), 7? = |(3-»), 7l=7 (1 _ 9) , 


and f k i, k = 0,1,2,3 is given in Eqs. (2.181 - (2.20) and (2.49), respectively. 


We deduce that 


(3.3) 


fi+i = 


fh | if0 = 1 ’ 


f 6 

0+5 


if 0 = 0. 


We then propose a new scheme in which / J+ i is chosen between /® +1 and /® +i adaptively. Hence, 
the scheme is 5th-order upwind or 6th-order central depending on the smoothness of the stencils S 5 and 
S 6 . We expect that this will get over the drawback of accuracy degeneration of the above mentioned 


central 6th-order schemes. To proceed, we first rewrite Eq. (3.1) using instead the non-linear weights 
w®’s as follows, 


(3.4) 

where, for k = 0,..., 3, and 

(3.5) 

Here, 8 


i — , io , , ,8 ll i,,Sf2 , , ,8 13 
fj+i - W o/ J + I +W 1 / 3 + i+14 2 / j+ i+U 3 / j+ i, 


Wfc = 




< = il\ i + 


£ + fik 


r” is the smoothness indicator of the large stencil, and ftk is the smoothness indicator of the 
sub-stencil Sk- We define these indicators in the following subsection. 

3.1. The Central Smoothness Indicators /?*. and New t . For a 6th-order central scheme over 
the large stencil S 6 , spurious oscillations are expected to occur near discontinuities. To overcome this, 
WENO-NW6 and WENO-CU6 choose to construct their (3 3 over all points of S 6 . A more careful 
observation reveals that this cost can be reduced in the following way. We remind the principle of 
WENO schemes is that there is at least one smoothest sub-stencil is used in the reconstruction of the 


numerical flux. We suppose that j3 3 follows Eq. (2.52), that is, it measures the smoothness of the 


most downwind S 3 only. We further assume that the grid size Aa; is so small that a discontinuity does 
not spread over two neighboring grid points, then for a 6th-order WENO scheme, the only case where 
oscillations occur is when a discontinuity is in between Xj and Xj+ 1 . In this case, both f°,i and f 3 ,i 


play main roles in the combination (3.41 since /3q and p 3 are much smaller than the other two. This leads 


to oscillations since the downwind f 3 i is wrongly chosen. To prevent this from happening, we choose 

2 _ 

P 3 to measure the smoothness of an extended sub-stencil S 3 := {xj,Xj+i,Xj+ 2 ,%j+ 3 } (see Fig. K^with 
S 3 extended by the dashed line). It is observed that S 2 is now a subset of S 3 and all sub-stencils share 
the point Xj. Hence, the case where oscillations occur is essentially eliminated. Moreover, the cost of 
computing p 3 is now much reduced since the computation involves only four grid points defined in S 3 
instead of six in S 6 as in WENO-NW6 or WENO-CU6. 

For 5th-order schemes, all sub-stencils are symmetric with respect to Xj. As a result, the smoothness 


indicators /3fc’s are also symmetric with respect to Xj (see Eqs. (2.291 - (2.31)). That is, (3q and P 2 are 
equal to each other up to order 0(Ax 4 ) in Taylor expansions. We recall that WENO discretizations 
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choose the sub-stencils depending on the non-linear weights uifs which are very sensitive to the smooth¬ 


ness indicators ftfs due to the latter’s smallness in smooth regions (see Eq. (2.271 for WENO-JS, Eq. 


(2.351 for WENO-Z, WENO-NW 6 , and WENO-CU 6 ). For the sensitivity, we mean that a small change 
in any 6 k leads to a large difference among afs, thus cofs. In that sense, the symmetry in terms of 
Taylor expansions of 6 k ’s reduces the effects of this sensitivity, especially in transition regions where the 
solution is smooth and discontinuous. We refer to Figs. [5] and [7] below for numerical evidences for this 
assessment in which the schemes with symmetric 6 fs (i.e., WENO-Z and WENO-0) show better results 
than the ones without this property. Unfortunately, the 6 th-order methods lack of this (comparing 63 


in Eq. (2.51) for WENO-NW 6 , and Eq. (2.55) for WENO-CU 6 with the other 6 fs, k = 0,1,2, defined 


Eqs. (2.29) - (2.31)). In our new scheme, we try to recover the property. We devise our new smoothness 
indicators in a central sense. That is, they are constructed based on the reconstruction polynomials 
which are symmetric with respect to 2 J+ i. In addition, it is shown below that the new indicators are 
symmetric in terms of Taylor expansions with respect to Xj. We notice that Taylor expansions about 


Xj are natural since the approximation of f{u) x is at the interval center Xj (see Eq. (2.3)) although the 


reconstruction of the numerical flux function h(x) is at the interface Xj + i (see Eqs. (2.6) and (2.44) for 
5th-order and 6 th-order schemes, respectively). 

Proceeding the reconstruction procedure as given in subsection 
with S 3 replacing S 3 , we obtain that 

2 / 7-2 - 7 /,_! + 11 fj , (fj-i - 3 /,•_! + 2 fj \ ^ , //,■_ 2 - 2/,-_! + fj 


but instead around 2 , , 1 =0 and 

J' 2 


(3.6) 


/>) = 


6 


Ax 


x + 


2Ax 2 


(3.7) 

(3.8) 

(3.9) 


f\x) = ^- 1+5/i + 2/j+1 + ( fj+ \ ~ fj | 2 


f\x) = 


V Ax 


2fj + 5/j+l — fj +2 , ffj+l-fj 


\ Ax 


fj-i ~ 2 f + fj+i \ 2 

2A2 2 1 


fj ~ 2 fj+l + fj +2 \ x 2 


2Ax 2 


f3( \ _ 3 fj + 13/7 + 1 - 5 /7+2 + fj +3 /-H/7 + 9/7 + 1 + 3/j+2 - fj+ 3 \ 

1 [X) ~ 12 + V 12Ax ) 

3fj - 7 /j+l + 5/7+2 - / 7+3 \ x 2 f - fj + 3/7 + 1 - 3/7+2 + fj+ 3 ^ ^3 


4A2 2 


6 A 2 3 


We notice that these reconstruction polynomials are different from those given in Eqs. (2.15) - (2.17) 
which are constructed symmetrically with respect to 2 j = 0. Substituting these into Eq. (2.28) with 
f k { 2 ) replacing f k { 2 ), k = 0,1,2,3, we deduce the new central smoothness indicators as follows, 


(3.10) 


(3.11) 


A) — Y^C/7-2 - 2/7-1 + fj ) 2 + { fj -2 - 3/7-1 + 2/7) 2 

= f 2 Ax 2 + ff" A2 3 + (|/" 2 - A2 4 + e>(A2 5 ), 

pi = ^(/i-l “ 2 /l + fj+i) 2 + (/7+I “ /j ) 2 

= f' 2 Ax 2 + f'f'Ax 3 + Q /" 2 + */ 7 "') A2 4 + 0 (A 2 5 ), 

62 = j^ifj ~ 2 fj+i + fj+2) 2 + (fj - fj+1) 2 

= f 2 Ax 2 + f'f'Ax 3 + Q /" 2 + i/ 7 "') A2 4 + 0 (A 2 5 ), 


(3.12) 
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and 


(3.13) 


13 


A — (3/j - 7/j+i + 5fj +2 - 7+3) + (2/j+i - 3 /j + 2 + /j+ 3 ) 

= f' 2 Ax 2 + f'f'Ax 3 + (\f" 2 - Iff") Ax 4 + 0( Ax 5 ), 


where the derivatives are evaluated at x = Xj. We note that for the most downwind /I 3 , we treated 

7(x) in Eq. (3.9) as a 2nd-degree polynomial by ignoring the 3rd-degree term when substituting it into 

2 


Eq. (2.28) so that (g|>/ 3 (a:)) 


dx is the highest-order variation. This is for the consistency with 

the other reconstruction polynomials f k ’ s, k = 0,1,2, which are of lower-degree. We also modified the 
second term in the obtained indicator so that its Taylor expansion agrees with that of /3q up to order 
0( Ax 4 ); thus all (3fs are now symmetric with respect to Xj in Taylor expansions, which is our goal in 
designing these new smoothness indicators. The original smoothness indicator (3$ for f3 3 obtained from 
Eq. (2.28) was as below, 

(3.14) fi 3 = —(3 f - 7f j+ i + 5/ j+ 2 - fj +3) 2 + + 9/j+i + 3/ i+2 - 7+ 3 ) 2 . 

In order to enhance the dispersion of the scheme, following the approach by Taylor et al. f |TWM07) 1. 
we set a restriction on the smoothness indicators as below, for k = 0,..., 3, 


(3.15) 

where 

(3.16) 


fik = 


0, if R(P) < a R , 
Pki otherwise; 


R(fi) = 


max fc (/3 fc ) 


£ + min k (0 k ) 

Here, a R is a threshold value depending on the configurations of flows. a R is taken small for flows with 
the presence of shocks. For a detailed discussion, consult [TWM07| . 

We next devise the smoothness indicator of the large stencil S e . Since the one proposed by Yamaleev 
and Carpenter in [ YC09 (see Eq. (2.53)) is too dispersive and may lead to oscillations (see the evidence 
in |YC091 1. we introduce a new smoothness indicator r e which is based on Eq. (2.28) but at a much 
higher order variations for the large stencil S 6 . We consider the following r’s, 


(3.17) 


and 


r 5 = Ax 5 


13 


r / 3 - 

^ /5(X) 


dx + Ax 7 




dx 


— ~ 4 fj-i + 6 fj - 4/j+i + fj+ 2 ) + (—fj -1 + 3 fj - 3/ J+ i + fj + 2 ) 


= f" ,2 Ax 6 + f'f^Ax 7 + Q/'"/ (5) + ^(/ (4) ) 2 ) Ax - 8 + (D(Ax 9 ), 


tq = Ax 7 


13, 


d 4 

s?™ 


dx + Ax 9 


d 5 ~ 


dx 


= ^(-7-2 + 57-1 - 10 f + I 0 f j+ 1 - 57+2 + fj +3 ) 

+ 4(7-2 - 37-1 + 2 /i + 2 /j'+i - 37+2 + fj+3) 2 
= (/ ( 4 ) ) 2 Ax 8 + fVf^Ax 9 + (J/( 4 >/ (5) + *(/< 5 >) 2 ) Ax 10 + ©(Ax 11 ), 


(3.18) 
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where f 5 (x) and / 6 (x) are the central reconstruction polynomials around x J+ i = 0 constructed in a 
similar way as with / fe (x)’s in Eqs. (3.6 1 - (3.9) but for the large stencils S 5 and S 6 , respectively; and 
the derivatives are evaluated at x,-. 


We then choose our t° and set 9 in Eq. (3.2) as follows 


(3.19) 


0 r e ,0) = 


(r 6 ,0) ifr 6 <r 5 , 
(r 5 ,l) ifr 6 >r 5 . 


It is noted that by choosing such t u and 9 as in Eq. (3.19), the scheme achieves a 6th-order in smooth 


regions since r 6 -C r 5 ; whereas it adaptively chooses the smoother large stencil between S 5 and S 6 
in the WENO reconstruction near discontinuities or unresolved regions. The new scheme now chooses 
the smoothest not only sub-stencils but also large one in the reconstruction procedure. The non-linear 


weights follow Eq. (3.5) above. Since our new method depends on 9 to switch between a 5th-order 


upwind and 6th-order central scheme, we name it WENO-0. In the numerical simulations below, we use 
the name WENO-06 for the compatibility with the other 6th-order schemes. 


Remark 3.1. 

i. Although the definition of r e has a switching mechanism by an if statement, it does not ruin 
the methodology of WENO schemes. This is because the switching applies to the smoothness 
indicator of the large stencil, not to the choice of smoother sub-stencils. 

ii. Although the switching is discontinuous in nature, WENO-0 is robust for problems with highly 
unstable fluid flows. We illustrate this by conducting a numerical simulation of the Rayleigh- 
Taylor instability problem in the below section. 

iii. The role of e is well investigated in [HAP05| . We note here that for cases with increasing number 
of vanishing derivatives, since both r and /?/. are very small at critical points, e does play roles 
to sustain the designed formal accuracy order. For this reason, except for WENO-JS where 
e = 10 -6 , we choose s = 10 -10 for other schemes. See the accuracy tests in the below section. 


In the next step, we test the accuracy, and resolutions of the new scheme. 


3.2. Accuracy Tests. We note that either 75 or tq is chosen in Eq. (3.19), the sufficient condition 
(2.50) is always satisfied. Hence the new scheme is 6 th-order in smooth regions. 

For the tests of accuracy, we choose the linear scalar conservation law, 


(3.20) 


u t + u x = 0, x G (-1,1), 
u(x, 0) = u 0 (x), 


subject to periodic boundary conditions. The following initial data are considered: 
• Initial condition 1: 


(3.21) 


uq(x) = sin(7rx); 


and 

• Initial condition 2: 
(3.22) 


uq(x) 



where fc — 1 is the number of vanishing spatial derivatives at x = — I, that is, 0 = ^ 
d^-Vf I / d w f I 

d x {h-l) la ,—0 / Q x ( fc ) x =fJ' 
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TABLE 1. Convergence of zt t + u x 


0 with initial conditions (3.21) and (3.22), at time t = 1. 




Eq. 1 

3.21 ) 

Eq. [3.22), 

k = 2 

Eq. ,'3.22), 

k = 3 


N 

L 1 error 

L°° error 

L 1 error 

L°° error 

L 1 error 

L°° error 

WENO- 

40 

4.5E-07 (-) 

3.4E-07 (-) 

4.5E-04 (-) 

2.2E-03 (-) 

4.2E-05 (-) 

1.3E-04 (-) 

CU6 

80 

6.9E-09 (6.0) 

5.4E-09 (6.0) 

3.8E-05 (3.6) 

2.0E-04 (3.5) 

8.9E-06 (2.2) 

4.5E-05 (1.5) 


160 

1.1E-10 (6.0) 

8.4E-11 (6.0) 

6.5E-07 (5.9) 

3.6E-06 (5.8) 

1.3E-07 (6.1) 

8.1E-07 (5.8) 


320 

4.1E-13 (8.0) 

3.8E-13 (7.8) 

1.1E-08 (5.9) 

6.0E-08 (5.9) 

1.8E-09 (6.1) 

1.1E-08 (6.2) 

WENO- 

40 

4.5E-07 (-) 

3.4E-07 (-) 

5.0E-04 (-) 

2.2E-03 (-) 

4.3E-05 (-) 

1.4E-04 (-) 

NW6 

80 

6.9E-09 (6.0) 

5.3E-09 (6.0) 

4.1E-05 (3.6) 

2.2E-04 (3.4) 

8.7E-06 (2.3) 

4.7E-05 (1.6) 


160 

1.1E-10 (6.0) 

8.4E-11 (6.0) 

6.4E-07 (6.0) 

3.6E-06 (5.9) 

1.4E-07 (5.9) 

9.6E-07 (5.6) 


320 

4.1E-13 (7.8) 

3.6E-13 (7.9) 

1.0E-08 (5.9) 

6.0E-08 (5.9) 

1.8E-09 (6.3) 

1.1E-08 (6.5) 

WENO- 

40 

4.5E-07 (-) 

3.4E-07 (-) 

3.7E-04 (-) 

1.8E-03 (-) 

3.2E-05 (-) 

1.2E-04 (-) 

(96 

80 

6.9E-09 (6.0) 

5.3E-09 (6.0) 

4.2E-05 (3.1) 

2.3E-04 (2.9) 

4.8E-06 (2.7) 

2.2E-05 (2.4) 


160 

1.1E-10 (6.0) 

8.4E-11 (6.0) 

7.6E-07 (5.8) 

3.9E-06 (5.9) 

1.2E-07 (5.3) 

6.3E-07 (5.1) 


320 

4.1E-13 (8.0) 

3.7E-13 (7.8) 

1.3E-08 (5.9) 

8.0E-08 (5.6) 

2.1E-09 (5.9) 

1.1E-08 (5.9) 


L 1 and L°° errors of 6th-order schemes at time t = 1 are measured and listed in Table [T] together 
with the order of accuracy (in brackets), and are plotted in Fig. [ 4 J We choose the time step At = Ax 6 / 3 
so that the numerical errors in time do not contribute to the results. In the figures, we also show the 
errors of the 5th-order schemes for comparison. It is observed that for all initial conditions, all schemes 
converge to the designed order of accuracy. Moreover, the errors of the new WENO-06 schemes are 
almost indistinguishable with those of other 6th-order ones, or even better at some grid sizes. 
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FIGURE 4. Convergence of Eq. (3.201 at time t = 1. Left: Initial condition (3.21); 
Middle: Initial condition (3.22) with k = 2. Right: Initial condition (3.22) with k = 3. 
Top: in L 1 norm; Bottom: in L°° norm. 


3.3. Resolution Tests. We now test if our new scheme overcomes the loss of accuracy of WENO-CU6 
and WENO-NW6. We revisit the initial condition given in example 0 which is as follows, 

(3.23) uq(x) = max(sin(7rx), 0), x £ (— 1,1). 


The numerical solution obtained from our new scheme is added and shown in Fig. [5j together with 
those given in example 0 We also plot the pointwise errors in L°° norm in the same figure. It is shown 
that WENO-06 approximates the critical region around x = —0.1 much better than WENO-CU6 and 
WENO-NW6. Indeed, the pointwise errors of the former around this region is comparable to those of 
WENO-JS and WENO-Z, where WENO-NW6 and WENO-CU6 show a loss of accuracy, which in turn 
causes problems for approximating solutions where symmetry is required. See below tests for numerical 
evidence. 

The nonlinear weights Wfc’s of these schemes are plotted in Fig. [6] We observe that around the critical 
point, the uq-’s of WENO-Z and WENO-06 are stable and converge to their optimal values q/ds. We 
note that for the latter scheme, the nonlinear weights keep fluctuating between the optimal weights of 
the 5th-order upwind and 6th-order central linear schemes. This clearly shows the effect of switching 
mechanism (3.19) in improving the accuracy of the scheme near the critical region. We also notify the 
non-convergence of Wfc’s of the WENO-CU6 and WENO-NW6 schemes around the critical region. This 
shows the improvement of our new scheme over the other 6tli-order ones. 
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FIGURE 5. Left: Numerical solutions of Eq. (2.1) with initial condition (3.23) at time 
t = 2.4. Grid 200. Middle: Zoom near the critical point. Right: Pointwise errors in log 
scale. 




WENO-NW6 



WENO-06 




FIGURE 6. Distribution of the non-linear weights for the initial data (3.23). From top 
to bottom, left to right: WENO-Z, WENO-CU6, WENO-NW6, and WENO-06. 


4. Numerical Results 

In this section, we perform a number of tests to compare the results of our new scheme with those 
obtained from the other WENO schemes, including the 5th-order upwind WENO-JS, WENO-Z, and the 
6th-order central WENO-CU6, WENO-NW6. 
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FIGURE 7. Left: Linear advection Eq. (3.20) with initial condition (4.1). Time t = 6.3. 
Grid 400. The others: zooms at critical regions. 


4.1. Scalar Conservation Laws. 


4.1.1. TEST 1: Linear Case. We solve the one-dimensional linear advection equation (3.20) with the 
following initial condition uq(x) which contains a C°° Gaussian, a square wave, a triangle, and a semi¬ 
ellipse (see | lJS961 1. 


(4.1) 


u 0 (x) = < 


5 [G(x, P,z-5) + 4 G{x, (3, z) + G(x, (3,z + < 5 )], 

1 , 

1 - |10(x — 0.1)|, 

^[F(x, a, a — S) + 4 F(x, a, a) + F(x, a,a + < 5 )], 

0 , 


-0.8 < x < -0.6, 
-0.4 < x < -0.2, 
0 < x < 0.2, 

0.4 < x < 0.6, 
otherwise, 


where 


(4.2) 


G{x, /3, z) = exp(-/3(z - zf), 


(4.3) F[x, a, a) = \J max(l — a 2 {x — a) 2 ,0); 

the constants are z = —0.7, 5 = 0.005, /? = a = 0.5, and a = 10. 

We compute the solution up to time t = 6.3 with N = 400 and periodic boundary conditions. The 
results obtained from the WENO-CU6, WENO-NW6, and WENO-06 schemes are plotted in Fig. [7| We 
choose an = 50. Zooms around the shocks and top of the semi-ellipse are also shown in the same figure. 
It is observed that WENO-06 is comparable to WENO-NW6 in capturing the shocks, but the former is 
much better than the latter and WENO-CU6 in approximating top of the semi-ellipse. 


4.1.2. TEST 2: Nonlinear Case. For this, we choose the Burgers equation 


(4.4) 


Ut 


= 0, x € (—1,1), 


u(x,0) = u 0 {x), 


subject to periodic boundary conditions. 

In Fig. [8j we show the numerical results of the 6th-order WENO schemes for the initial condition 


u 0 (x) = - 


(4.5) 


sin(7ra:); 
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FIGURE 8. Burgers’ Eq. (4.4). Grid 200. Left: initial condition (4.5) at time t = 1.5; 
Right: initial condition (4.6) at time t = 0.55. 


at time t = 1.5 and 

(4.6) uq(x) = - + sin( 7 ra:) 

at t = 0.55. We choose a grid of N = 200 grid intervals and an = 50. It is shown that the shocks are 
very well captured by all schemes. 

4.2. Euler Equations of Gas Dynamics. In this subsection, we consider the one-dimensional Euler 
equations of gas dynamics given in the following, 

(4.7) u t +f(u) x = 0, 
where 

(4.8) u = (p,pu,E) T , f( u ) = (pu,p + pu 2 , (E + p)u) T , 

where p , u, p , E are density, velocity, pressure, and total energy, respectively. The equation of state is 
as follows, 

(4.9) p=('y-l)^E- \pu 2 ^j , 7 = 1-4, 

where 7 is the ratio of specific heats. More details on the Euler equations can be found in, e.g., ILV921 . 

ITo97| . 

For all below numerical simulations, we apply the WENO schemes in characteristic fields of the flux 
f(u). That is, we first find an average of the Jacobian A^i at the interface x = Xj + 1 . For this, we apply 
the Roe’s mean matrix (see [RoSin . Then the eigenvalues A s ’s, L = [lsj^Li, R = [ r s]?Li the complete 
sets of the left and right eigenvectors, respectively, of Aj + 1 are determined. We next project the flux f(u) 
into the characteristic fields by left multiplying it with L. WENO schemes with a global Lax-Friedriclis 
flux splitting are applied to approximate the components of the flux. After that, the approximation in 
each characteristic field is projected back to the component space by a right multiplying with the matrix 
R. 
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FIGURE 9. Left: Sod’s problem with initial data (4.10). Time t = 1.7. Grid 300. 
Right: zoom at the contact discontinuity. 


4.2.1. TEST 3: Riemann Problems. We consider the shock-tube problems which are Eq. (4.7) with 
Riemann initial data. In particular, the Sod problem, the Lax problem, and the 123 problem are given 
below. 

• Sod’s problem: 


(4.10) ( P,u,p) = 

and the final time t = 1.7. 

• Lax’s problem: 


(0.125, 0, 0.1), —5 < x < 0, 

(1, 0, 1), 0 < a; < 5; 


(4.11) ( P,u,p) = 

and the final time t = 1.3. 

• 123 problem: 


(0.445, 0.698, 3.528), -5 < x < 0, 

(0.5, 0, 0.571), 0 < x < 5; 


(4.12) 


(p,u,p) = 


(1, -2, 0.4), —5 < £ < 0, 
(1, 2, 0.4), 0 < x < 5; 


and the final time t = 1. 

We apply a transmissive condition at both boundaries. The exact solution of these shock-tube prob¬ 
lems can be found in, for example, [Tb97l . For Lax’s and the blast waves problems, we choose an = 10; 
and an = 50 for other problems. 

Numerical results of the density obtained from all WENO schemes with a grid of ./V = 300 are shown 
in Figs. 0-0 respectively. We observe that for Sod’s and Lax’s problems, there are overshoots at 
the contact discontinuities for WENO-CU6, whereas WENO-06 gives the sharpest capturing without 
generating oscillations. For the 123 problem, WENO-CU6 shows the most wiggling behavior around the 
trivial contact discontinuity. 



























24 


WEIGHTED ESSENTIALLY NON-OSCILLATORY SCHEMES 




FIGURE 10. Left: Lax’s problem with initial data (4.11). Time t = 1.3. Grid 300. 


Right: zoom at the contact discontinuity. 




FIGURE 11. Left: The 123 problem with initial data (4.12). Time t = 1.0. Grid 300. 
Right: zoom at the trivial contact discontinuity. 


4.2.2. TEST Shock Density Wave Interaction, Shu-Osher’s test. We consider the following initial 
data, 


(4.13) 


(p,u,p) = 


(3.857143, 2.629369, 31/3), -5 < x < -4, 

(1 + 0.2sin(5x), 0, 1), -4 < x < 5, 


with zero-gradient boundary conditions. 

The problem simulates the interaction of a right-moving Mach 3 shock with a wavelike perturbed 
density whose magnitude is much smaller than the shock. As a result, a flow field of compressed and 
amplified wave trails is created right behind the shock. For more details, see [JS96 . In Fig. 12 we show 
the numerical results of the 6th-order WENO schemes at time t = 1.8 with grids of N = 200 and N = 400 
intervals. The “exact” solution is computed by WENO-JS with a fine grid N = 4000. It is shown that 
all schemes give satisfactory approximations of the compressed wavelike structures behind the shock. A 
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FIGURE 12. Shu-Osher’s problem with initial data (4.13). Density. Time t = 1.8. 
Left: medium grid 200; Right: fine grid 400. 


careful observation reveals that WENO-06 resolves the wave package as well as WENO-CU 6 , whereas 
WENO-NW 6 is more dissipative for both grid levels. 


4.2.3. TEST 5: Two Interacting Blast Waves. In this test, we show that our new scheme WENO-06 
passes the tough test of two interacting blast waves which the initial data are given as follows, 



f (1, 

0 , 

1000 ), 

0 < x < 0 . 1 , 

(4.14) 

e 

SL 

II 

0 , 

0 . 01 ), 

0.1 < X < 0.9 


{ (i, 

0 , 

100 ), 

0.9 < x < 1, 


and a reflective condition is applied at both boundaries. This problem is used to test the robustness of 
shock-capturing methods since many interactions are observed in a small area. A detailed discussion of 
this problem can be found in |WC84| . 

Numerical results of 6 th-order WENO schemes are computed up to time t = 0.038 with a grid of 
N = 801 and plotted in Fig. [13] for the density. The exact solution is approximated by WENO-JS with 
a much fine grid N = 4001. It is shown that all schemes well capture the shocks as well as contact 
discontinuities. A zoom near x = 0.745 indicates that WENO-06 gives better resolution than WENO- 
NW 6 and WENO-CU 6 . We also emphasize that there exists a stair-casing phenomenon in the solutions 
of the latter methods in this region, which is similar to that at the top of the semi-ellipse in test 1 (see 
Fig. [7]), and the 123 problem (see Fig. 11). 


4.3. Two-dimensional Euler’s Equations. In this subsection, we extend the problem to two-dimensional 
cases. We choose the 2D Euler equations which are as follows, 


(4.15) 


U t +f(u) K + g(u) y = 0, 
u(x,y,0) = u 0 (x,y), 


where u = (p, pu, pv, E) T , f(u) = (pu,p + pu 2 , puv,u(E + p)) T , g(u) = (pv, puv,p + pv 2 ,v(E + p)) T . 
The relation of pressure and conservative quantities is through the equation of state 


(4.16) 


P = (7 - 1 ) [e 



Here, we choose the ratio of specific heats 7 = 1.4. 
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FIGURE 13. Two interacting blast waves with 
t = 0.038. Grid 801. 


initial data (4.14). Density. Time 


4.3.1. TEST 6: Rayleigh-Taylor Instability. In the following tests, we show numerical evidence that 
WENO-06 maintains symmetry in the solutions much better than the other 6th-order schemes, and 
outperforms 5th-order schemes in resolving small-scaled structures occurring in flow configurations. We 
first simulate the Rayleigh-Taylor instability. The instability occurs where there is a heavy fluid falling 
into a light fluid (see {GGL088] , Ktj . [FSTY14j ). Following 2t], we set up the problem as follows. 
The domain is (x,y) = (—0.25,0.25) x (—0.75,0.75). Initial density has a discontinuity at the interface, 
i.e., p = 2 for y > 0 and p = 1 for y < 0. The pressure is set at hydrostatic equilibrium initially 
p = 2.5 — pgy where g = 0.1 is the gravitational acceleration. The ^-component velocity u = 0, while 
the y-component is perturbed with v = ^p(l + cos(47ro;))(l + cos(|7ry)) for a single mode perturbation. 
Boundary conditions are set periodic in ^-direction, and reflective in y-direction. The ratio of specific 
heats 7 = 1.4. We add — gp and — gpv in the y-momentum and energy equations of (4.15) as source 
terms. 

In Fig. [l4j we plot the density with 20 equally spaced contours obtained from 5th- and 6th-order 
WENO schemes at time t = 9.5 with a 120 x 360 grid. It is shown that the 6th-order schemes have much 
better numerical resolution comparing with the 5th-order ones. We notice that WENO-06 preserves 
the symmetry of the solution; whereas WENO-NW6 and WENO-CU6 do not. We conjecture the lack 
of symmetry of WENO-NW6 is due to the loss of accuracy around critical regions which is shown in 
previous numerical tests. The test also shows that the discontinuous switching of t 6 in Eq. (3.19) does 
not affect the robustness of the new WENO-0 scheme, even for problem with highly unstable fluid flows 
as the Rayleigh-Taylor instability. 
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WENO-Z: density WENO-NW6: density WENO-CU6: density WENO-06: density 






FIGURE 14. The Rayleigh-Taylor instability. Density at time t = 9.5. Grid 120 x 360. 
CFL = 0.5. From left to right: WENO-Z, WENO-NW6, WENO-CU6, and WENO-06. 


4.3.2. TEST 7: Implosion problem. The next numerical test is the implosion problem (see [At] . [ LW04] ) 
with initial data as follows, 


(4.17) 



(1,1) for x + y>\, 
(0.125,0.14) otherwise, 


and zero velocity everywhere initially. We choose reflecting conditions for all boundaries. 

Symmetry is important for this test. For such a scheme, due to the interactions of shock waves and 
reflecting boundaries, jets along the diagonal are created. Longer and narrower jets are produced for less 
dissipative schemes. 

In Fig. 15 we show the results obtained from different schemes on computational domain (0,1) x (0,1) 
at final time t = 5. We choose a grid of 400 x 400. For WENO-0, we choose an = 1- It is shown that 
only WENO-Z and WENO-0 well preserve the symmetry of the problem; whereas the other 6th-order 
schemes do not. The jets created by WENO-NW6 and WENO-CU6 tend to diverge from the main 
diagonal x = y. We also note that the jets produced by WENO-0 is much longer and narrower than 
those of WENO-Z, which means that the former scheme is less dissipative than the latter one. 


4.3.3. TEST 8: 2D Riemann Initial Data. The 2D Riemann problem is set up by assigning different 
constant states of (pfc,Mfc, Vk,Pk), k = 1,2,3,4, to four quadrants of the computational domain O = 
(0,1) x (0,1). The constant states are chosen so that there is only a single elementary wave, namely, 
shock-, rarefaction-, and contact-wave, connecting two neighboring quadrants (see [SGG93j b For our 
test, we choose the following configuration for the initial data, respectively, for quadrants 1, 2, 3, 4, 


(0.5313, 0, 0, 0.4), 
(1, 0.7276, 0, 1), 
( 0 . 8 , 0 , 0 , 1 ), 
k (1, 0, 0.7276, 1), 


x > 0.5, y > 0.5, 
x < 0.5, y > 0.5, 
x < 0.5, y < 0.5, 
x > 0.5, y < 0.5, 


(4.18) 


(p,u,v,p) = < 
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WENO-Z: density 



0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 

WENO-CU6: density 




FIGURE 15. The implosion problem. Density with 20 contours uniformly distributing 
from 0 to 1. Grid 400 x 400. Final time t = 5. Left to right, top to bottom: WENO-Z, 
WENO-NW6, WENO-CU6, WENO-06. 


which has shocks through quadrants 1-2 and 1-4, and contact discontinuities through quadrants 2-3 
and 3-4. Transmissive boundary conditions are imposed on all boundaries for these two cases. 


The approximations of the density with initial data (4.18) at time t = 0.25 are plotted in Fig. 16 


with 50 contours for WENO-NW6 and WENO-06. Here, we use a fine grid with 1000 x 1000 intervals 
for the capturing of the small vortices along the contacts. Zooms near the spirals region are also shown 
on the right column in the same figure. Again, we observe a better performance of the WENO-06 over 
WENO-NW6 and WENO-CU6 schemes over these small-scaled structures without oscillations on the 
contours. 


4.3.4. TEST 9: Double Mach Reflection of a Strong Shock. Finally, we investigate the double Mach 
reflection of a strong shock which is a typical benchmark test for shock-capturing methods. The problem 
simulates the reflection occurring when a simple planar shock interacts with a wedge making with the 
x-axis an angle a. The strength of the moving shock is characterized by the Mach number M s . For a 
double Mach reflection problem, M s = 10 and the wedge angle is chosen as a = 30°. Detailed discussions 
on this type of problems can be found in IWC84j and the references therein. For numerical purpose, we 
choose the computational domain H = (0,4) x (0,1). Initially the shock is located at xq = 1/6, inclined 
with the 2 -axis by the angle 90° — a. Inflow and zero gradients conditions are imposed on the left and 
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WENO-NW6: density 


WENO-NW6: density 



FIGURE 16. Left: The 2D Riemann problem with initial data (4.18). Density with 50 
contours. Time t = 0.25. Grid 1000 x 1000. Right: Zoom at the spirals. From top to 
bottom, respectively: WENO-NW6, WENO-CU6, WENO-6>6. 


right boundaries, respectively. On the bottom one, a reflective condition is applied to the interval [xo,4] 
representing the wedge, and the exact post-shock state is imposed over [0, xo]. The top boundary is 
treated in a way that there are no interactions of the shock with this boundary. That is, the exact 
post- and pre-shock states are employed over the intervals [0, x s (i)] and [x s (i),4], respectively, on the 


top boundary. Here, x s (t) = Xq + 


1 


A/. o n 


tan 60° cos 30° 


t, where a pre is the sound speed of the pre-shock 


























30 


WEIGHTED ESSENTIALLY NON-OSCILLATORY SCHEMES 


state, is the location of the shock in time. These states can be computed exactly when one of them is 
pre-described (see, e.g., [To97] 1. In particular, for our problem the initial conditions are given as follows, 


(4.19) 


(p,u,v,p) = 


( 8 , 8.25cos30°, — 8.25sin30°, 116.5), x < xq + tan y 60o , 
(1.4, 0, 0, 1), x > x 0 + tan y 6QO . 


Numerical results of the density obtained from the 5tlr-order WENO-Z, the 6 th-order WENO-NW 6 , 
WENO-CU 6 , and WENO-06 schemes at time t = 0.2 are plotted in Fig. [IT] with 30 contours. For this 
case, we choose a fine grid of 800 x 200 points for all schemes. We notice the rendering of small vortices 
at the end of the slip line and the wall jet, starting from WENO-Z and becoming clearer for the 6 th-order 
schemes. The zoom-in on the Mach stems region shown in Fig. [18] reveals that the WENO-06 scheme 
gives more satisfactory resolution than the WENO-NW 6 and WENO-CU 6 ones. 


5. Conclusion 

In this work, we have presented a new WENO-0 scheme which adaptively switches between a 5th-order 
upwind and 6 th-order central scheme, depending on the smoothness of not only the sub-stencils but also 
the large one. Unlike the other 6 th-order WENO methods in which this switch depends solely on the 
smoothness of the most downwind sub-stencil, it is that of the large stencil which decides this mechanism 
in our new scheme. Main features of our new scheme are that the new scheme combines good properties 
of both 5th-order upwind and 6th-order central schemes. That is, the new scheme is more dispersive than 
the 5th-order ones in terms of better resolution of small-scaled structures and capturing discontinuities. 
Moreover, the scheme overcomes the loss of accuracy around some critical regions and has the ability to 
maintain symmetry in the solutions which are drawbacks of other comparing 6th-order WENO schemes. 

We have also developed the new smoothness indicators of the sub-stencils 0 /s which are symmetric 
in terms of Taylor expansions around the point Xj and a new r 0 for the large stencil. The latter is 
chosen as the smoother one among two candidates which are computed based on the possible highest- 
order variations of the reconstruction polynomials in L 2 sense. From then, value of the parameter 9 is 
determined to decide if the scheme is 5th-order upwind or 6 th-order central. 

A number of numerical tests for both scalar cases, linear and nonlinear, and system case with the 
Euler equations of gas dynamics are carried out to check the accuracy, resolution, and robustness of 
our new scheme. It is shown that our new method is more accurate than WENO-JS, WENO-Z, and 
WENO-CU 6 ; more robust than WENO-CU 6 and WENO-NW 6 ; and outperforms comparing schemes in 
capturing small-scaled structures, and around critical regions. 

Numerical simulations of higher dimensional problems will be investigated in a subsequent work. 
Since the new smoothness indicators 0 3 and r® are constructed in a systematic manner, we expect a 
development of the scheme to a higher order of accuracy. This will be considered in our future work. 
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WENO-Z: density 



WENO-NW6: density 



Figure 17. 


The double-Mach reflection problem with initial data (4.191. Density with 


30 contours. Time t = 0.2. Grid 800x200. From top to bottom, respectively: WENO-Z, 


WENO-NW6, WENO-CU6, and WENO-06. 
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WENO-Z: density 


WENO-NW6: density 



Figure 18. 


The double-Mach reflection problem with initial data (4.19). Zoom at the 


double Mach stems region. From top to bottom, respectively: WENO-Z, WENO-NW6, 


WENO-CU6, and WENO-6>6. 
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